### interesting spp
library(raster)
## Loading required package: sp
source('https://raw.githubusercontent.com/oharac/src/master/R/common.R')
## ── Attaching packages ─────────────────────────────────────── tidyverse 1.3.0 ──
## ✓ ggplot2 3.3.3 ✓ purrr 0.3.4
## ✓ tibble 3.1.0 ✓ dplyr 1.0.5
## ✓ tidyr 1.1.3 ✓ stringr 1.4.0
## ✓ readr 1.4.0 ✓ forcats 0.5.1
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## x tidyr::extract() masks raster::extract()
## x dplyr::filter() masks stats::filter()
## x dplyr::lag() masks stats::lag()
## x dplyr::select() masks raster::select()
## here() starts at /home/ohara/github/bd_chi
source(here::here('common_fxns.R'))
library(sf)
## Linking to GEOS 3.8.0, GDAL 3.0.4, PROJ 7.0.0
knitr::opts_chunk$set(fig.width = 6, fig.height = 4, fig.path = 'figs/',
echo = TRUE, message = FALSE, warning = FALSE)
# cellid_mol <- raster('_spatial/cell_id_mol.tif')
# ocean_mol <- raster('_spatial/ocean_area_mol.tif')
# values(cellid_mol)[is.na(values(ocean_mol))] <- NA
# blank_latlong <- raster(ext = extent(c(-180, 180, -90, 90)), crs = '+init=epsg:4326', res = .1)
# cellid_latlong <- projectRaster(cellid_mol, blank_latlong, method = 'ngb')
# writeRaster(cellid_latlong, '_spatial/cellid_latlong.tif', overwrite = TRUE)
cellid_latlong <- raster('_spatial/cellid_latlong.tif')
ca_cellid_latlong <- crop(cellid_latlong, extent(c(-130, -115, 30, 42)))
ca_cells <- values(ca_cellid_latlong) %>% .[!is.na(.)]
spp_ids <- get_incl_spp()$iucn_sid %>% unique()
spp_cells_list <- parallel::mclapply(spp_ids, mc.cores = 40, FUN = function(id) {
f <- sprintf(file.path(dir_bd_anx, 'spp_rasts_mol_2020/iucn_sid_%s.csv'), id)
x <- data.table::fread(f)
})
spp_cells_df <- spp_cells_list %>%
setNames(spp_ids) %>%
bind_rows(.id = 'iucn_sid') %>%
mutate(iucn_sid = as.integer(iucn_sid))
ca_spp_cells <- spp_cells_df %>%
filter(cell_id %in% ca_cells)
ca_spp <- ca_spp_cells %>%
select(iucn_sid) %>%
distinct() %>%
left_join(get_incl_spp())
n_ca_spp <- ca_spp %>%
select(iucn_sid, sciname, comname, cat_score, assess_gp) %>%
distinct() %>%
group_by(assess_gp) %>%
mutate(n_gp = n()) %>%
left_join(get_spp_range() %>% filter(eez == 0))
### A few identified species
spp_of_note <- c(46967817, 46967863, 8005, 39385, 19488,
17028, 4162, 8099, 22697742,
21858 , 21860, 170341, 193734, 10088)
imps <- get_incl_spp() %>%
filter(iucn_sid %in% spp_of_note) %>%
select(iucn_sid, sciname, comname, cat_score) %>% #, assess_gp, desc, taxon) %>%
distinct() %>%
left_join(read_csv('_output/imp_range_by_spp_2013.csv'), by = 'iucn_sid') %>%
filter(impact_km2 > 0) %>%
mutate(imp_pct = round(impact_km2 / range_km2 * 100, 3),
inc2_pct = round(incr2_km2 / range_km2 * 100, 3),
inc3_pct = round(incr3_km2 / range_km2 * 100, 3),
dec2_pct = round(decr2_km2 / range_km2 * 100, 3),
dec3_pct = round(decr3_km2 / range_km2 * 100, 3)) %>%
select(-ends_with('km2'), range_km2) %>%
distinct()
Focus on turtles and whale sharks
spp_df <- get_incl_spp() %>%
filter(iucn_sid == 19488 | assess_gp == 'sea_turtles') %>%
select(-code) %>%
distinct()
results_df <- read_csv('_output/imp_range_by_spp_2013.csv') %>%
filter(iucn_sid %in% spp_df$iucn_sid) %>%
mutate(imp_pct = round(impact_km2 / range_km2 * 100, 3),
imp_2plus_pct = round(impact_2plus_km2 / range_km2 * 100, 3),
incr2_pct = round(incr2_km2 / range_km2 * 100, 3),
incr3_pct = round(incr3_km2 / range_km2 * 100, 3),
decr2_pct = round(decr2_km2 / range_km2 * 100, 3),
decr2_pct = round(decr3_km2 / range_km2 * 100, 3)) %>%
filter(impact_km2 > 0) %>%
select(-ends_with('km2'), range_km2) %>%
distinct() %>%
left_join(spp_df %>%
select(-stressor, -category) %>% distinct())
DT::datatable(results_df)
cell_id_rast <- raster('_spatial/cell_id_mol.tif')
land_sf <- read_sf('_spatial/ne_10m_land/ne_10m_land.shp') %>%
st_transform(crs(cell_id_rast))
for(spp in unique(spp_df$iucn_sid)) {
### spp <- spp_df$iucn_sid[1]
### spp <- 19488
tmp_df <- spp_df %>%
filter(iucn_sid == spp)
spp_name <- tmp_df$sciname[1]
spp_name <- str_remove(spp_name, 'Ocean|ulation')
spp_name2 <- tmp_df$comname[1]
message('Processing ', spp_name)
spp_range_f <- sprintf(file.path(dir_bd_anx, 'spp_rasts_mol_2020', 'iucn_sid_%s.csv'), spp)
cells <- read_csv(spp_range_f)
range_rast <- raster(cell_id_rast)
range_rast[values(cell_id_rast) %in% cells$cell_id] <- 0
strs <- tmp_df$stressor %>% unique()
str_df <- lapply(strs, FUN = function(str) {
### str <- strs[1]
str_f <- sprintf(file.path(dir_bd_anx, 'spp_str_rasts/%s/spp_intsx_%s_%s.csv'),
str, str, spp)
df <- read_csv(str_f) %>%
filter(year == 2013) %>%
mutate(stressor = str,
impact = 1)
}) %>%
bind_rows()
for(str in strs) {
### str <- strs[1]
df <- str_df %>%
filter(stressor == str) %>%
as.data.frame()
r <- range_rast
r[values(cell_id_rast) %in% df$cell_id] <- 1
r_df <- rasterToPoints(r) %>%
as.data.frame() %>%
setNames(c('x', 'y', 'val'))
p <- ggplot() +
ggtheme_plot() +
theme(axis.title = element_blank()) +
geom_sf(data = land_sf, fill = 'grey70', color = 'grey30', size = .1) +
geom_raster(data = r_df, aes(x, y, fill = val)) +
scale_fill_viridis_c() +
labs(title = sprintf('%s %s %s', spp_name2, spp_name, str))
print(p)
}
sum_df <- str_df %>%
group_by(cell_id) %>%
summarize(n_strs = n())
sum_rast <- raster::subs(cell_id_rast, sum_df, by = 'cell_id', which = 'n_strs')
values(sum_rast)[is.na(values(sum_rast))] <- values(range_rast)[is.na(values(sum_rast))]
sum_df <- rasterToPoints(sum_rast) %>%
as.data.frame() %>%
setNames(c('x', 'y', 'val'))
p <- ggplot() +
ggtheme_plot() +
theme(axis.title = element_blank()) +
geom_sf(data = land_sf, fill = 'grey70', color = 'grey30', size = .1) +
geom_raster(data = sum_df, aes(x, y, fill = val)) +
scale_fill_viridis_c() +
labs(title = sprintf('%s %s cum', spp_name2, spp_name))
print(p)
}
